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Abstract 

The Yilgarn Banded Iron Formations of Western Australia are topographical features that behave as terrestrial islands within 
the otherwise flat, semi-arid landscape. The formations are characterised by a high number of endemic species, some of 
which are distributed across multiple formations without inhabiting the intervening landscape. These species provide an 
ideal context for phylogeographic analysis, to investigate patterns of genetic variation at both spatial and temporal scales. 
We examined genetic variation in the spirostreptid millipede, Atelomastix bamfordi, found on five of these Banded Iron 
Formations at two mitochondrial loci and 11 microsatellite loci. Strong phylogeographic structuring indicated the five 
populations became isolated during the Pleistocene, a period of intensifying aridity in this landscape, when it appears 
populations have been restricted to pockets of moist habitat provided by the formations. The pattern of reciprocal 
monophyly identified within the mtDNA and strong differentiation within the nuclear microsatellite data highlight the 
evolutionary significance of these divergent populations and we suggest the degree of differentiation warrants designation 
of each as a conservation unit. 
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Introduction 

Islands provide a unique context for studying the evolutionary 
processes that drive species distributions. They are natural 
laboratories, allowing examination of processes of speciation, 
and patterns of colonisation, persistence, dispersal and extinction 
[1]. The genetic consequences of island distributions are well 
documented by phylogeographic and population genetic studies 
that infer historical patterns of natural and sexual selection, drift 
and their interactions with historical and ongoing dispersal [2—5]. 
The study of terrestrial or 'habitat' islands appears to offer a 
different challenge to oceanic islands, as the matrix between 
islands can be more permeable than water barriers, resulting in 
potentially complex spatial and temporal patterns of genetic 
structure [6,7]. The most conspicuous and well-studied terrestrial 
islands include mountain Ranges that form 'sky islands', such as 
the tepui summits of South America, the Western Ghats of 
southern India and the Madrean and Rocky Mountain sky island 
systems of central and northern America [6,8-12]. Here, 'island 
species' occupy disjunct populations on either mountain tops or 
the lowland valleys separating mountains. Other, less conspicuous 
topographic habitat islands are found throughout the world and 
include granite outcrops, inselbergs and Banded Iron Formations 
(BIFs), many of which display high endemism. Although 
topographically less dramatic, these provide equally compelling 
bases for phylogeographic study [13-17]. 



Inferring the evolutionary history of terrestrial island endemics 
can be difficult given that isolation is not always absolute, either 
spatially and temporally. Dispersal events can occur across valleys 
and lowlands for mountain-top endemics or over mountains for 
valley floor endemics, depending on the vagility of the organism 
and their environmental requirements [6,7]. Despite this, 
phylogeographic patterns, although potentially complex, may 
reveal the degree of connectivity or isolation experienced by island 
populations. 

Restricted gene flow between island populations over long 
periods of time, will result in substantial genetic divergence due to 
genetic drift and/or, local adaptation [18]. Given sufficient time, 
populations may show patterns of reciprocal monophyly, where all 
individuals can be traced to a common ancestor within that 
population [19,20]. Recurrent migration and gene flow between 
islands on the other hand, has a homogenising effect on genetic 
structure unless there is strong selection favouring local genotypes 
on each island. For example, populations of yellow bellied 
marmots, restricted to the Great Basin mountain-tops of North 
America, were believed to have been marooned on the rocky 
outcrops of mountain tops following habitat contraction during the 
Pleistocene [21]. A pattern of isolation by distance and high 
genetic diversity showed that occasional dispersal events across the 
arid scrublands of the Basin have been an important process 
connecting mountain-top populations [21]. A contrasting pattern 



PLOS ONE | www.plosone.org 



1 



March 2014 | Volume 9 | Issue 3 | e93038 



Phylogeography of a Semi-Arid Millipede 



was found in the tree Eucalyptus caesia, endemic to granite outcrop 
'islands' of Western Australia. This species has high levels of 
genetic divergence with little geographic pattern to the distribution 
of chloroplast haplotypes, reflecting long-term isolation and 
persistence without connectivity and emphasising the role of 
genetic drift in this landscape [22]. 

The Yilgarn Banded Iron Formations (BIFs) of south-western 
Australia are ancient, island-like rock outcrops that provide an 
ideal system in which to investigate island population dynamics. 
They border the transitional rainfall zone, a biogeographic region 
separating the high rainfall forest region from the arid interior and 
are characterised by high species endemism and richness [14,23]. 
BIFs vary widely in height (highest point approx. 700 m above sea 
level), shape and distance from one another and are comprised of 
banded iron talus slopes and areas of weathered duricrust [24] . 
These topographically complex formations occur in an otherwise 
flat landscape that is dominated by clay, silt and sand plains. The 
region has been geologically stable, un-glaciated and above sea 
level since the Permian (250 mya) [25-27]. Nonetheless there have 
been substantial changes in climate and vegetation communities 
over this time. At the beginning of the Genozoic (65 mya) the 
region was warm, humid and dominated by temperate rainforest 
[28]. Over time, the landscape became increasingly arid, 
beginning approximately 1 1 mya in the mid-Miocene, shown by 
the absence of regular flows in paleo-drainage systems [28]. 
Another aridification event occurred 7 mya in the late Miocene 
with several more throughout the Pliocene, shown by changes in 
vegetation structure in palynological deposits to more arid-adapted 
forms [28,29]. By the early Pleistocene (approx. 2.6 mya), the 
modern climate had more or less become established but 
fluctuations throughout this period continued, associated with 
successive glacial and interglacial periods, the latter with 
intensifying arid conditions [28]. 

Historically, the altitudinal relief of BIFs and their complex 
structure, including deep valleys and shaded gullies as well as 
exposed outcrops and talus slopes [24], would have provided an 
abundance of niche habitats that allowed for the persistence of 
both mesic and xeric-adapted species during changing climates 
[30]. Whilst this may account for the high species diversity found 
on the BIFs today, the high endemism [14,23,31] also suggests 
long periods of isolation, allowing for speciation and the 
development and persistence of a unique flora and fauna. The 
Yilgarn BIF endemics consist of species that are localised (found on 
one BIF), or regional (located on several BIFs but not in the 
intervening landscape) [23]. Regional endemics provide an 
opportunity to study phylogeographic patterns, to determine what 
processes have caused current distributions, and historically, how 
diversity has evolved and accumulated. The floral assemblages of 
the Yilgarn BIFs are well known, with many examples of regional 
endemics [23]. Although less is known of the faunal assemblages, 
taxa that exhibit short-range endemism consistent with limited 
dispersal capacity, such as millipedes, mygalomorph spiders, land 
snails and troglofauna [32-34], may also exhibit patterns of 
regional endemism, persisting in micro-habitats provided by the 
BIF. 

We investigated phylogeographic pattern in the spirostreptid 
millipede Atelomastix bamfordi, known to occur on five BIFs on the 
Yilgarn Plateau in Western Australia. Millipedes lack the waxy 
cuticle that protects other arthropods from desiccation and are 
reliant on moist environments for survival [32,35]. Atelomastix from 
the family Iulomorphidae, is a predominantly south-western 
Australian genus with many species distributed across the wetter 
parts of the southern and south-western coastline and mesic inland 
Ranges, including the Stirling and Porongorup Ranges [36,37]. 



Almost all species have restricted distributions and have been 
characterised as short-range endemics (geographic range 
< 10000 km 2 ) [32]. Atelomastix bamfordi is the most easterly 
occurring member of this genus. It inhabits moist habitats 
provided by BIFs, such as the base of rocky slopes, under large 
rocks and under dense leaf litter. The life cycle is probably annual 
with adult millipedes dying prior to, or during the harsh summer 
months (pers. comm. M Harvey). The species has been reported as 
abundant on the large Koolyanobbing Range [33] but appears 
scarcer on the drier BIFs to the north and east. Increased 
exploration and mining of iron ore within the region [23] has 
resulted in the need for a delicate balance between mining 
activities and conservation of BIF endemic species. Understanding 
the patterns and evolution of diversity in these sites may guide two 
key conservation issues: i) the viability and conservation of highly 
endemic species ii) the conservation of landscape features which 
maximise the chances of maintaining patterns of dispersal, gene 
flow and local adaptation that have shaped the history of extant 
diversity and continue to be major factors in the persistence of this 
diversity. 

Although unconnected today, BIF populations of A. bamfordi 
may have repeatedly come into contact through population 
expansion associated with historical climate change in this 
landscape. If this were the case, we would expect periods of gene 
exchange between neighbouring BIF populations to result in a 
pattern of isolation-by-distance, where proximate populations are 
more genetically related, or, if gene exchange has been extensive, 
to genetic homogeneity of populations [38]. Alternately, the 
current distribution may be due to contraction of a formerly 
widespread ancestor followed by long-term isolation and persis- 
tence on BIF islands. If this were the case, we would expect to see 
signatures of deep genetic divergence amongst BIF populations but 
no clear pattern of isolation by distance. Haplotypes would be 
conserved and there may be evidence of reciprocal monophyly, if 
sufficient time for lineage sorting has occurred [39]. This is the first 
phylogeographic study to examine patterns of population level 
divergence in a regional BIF endemic species from this important, 
biodiverse region. 

We evaluated the genetic variation present in A. bamfordi to 
investigate evolutionary relationships and patterns of genetic 
diversity and structure across BIF populations. Our goals were: 

1 . To determine patterns of genetic differentiation among the five 
known BIF populations and to use those data to infer the 
impact of historical processes on the distribution of millipedes 
across the landscape. 

2. To assess genetic diversity and structure to determine how 
limited dispersal capabilities affects structure within BIF 
millipede populations. 

3. To use patterns of genetic differentiation, both within and 
among BIFs to inform conservation and management of this 
species. 

Materials and Methods 

Sampling and DNA extraction 

A. bamfordi had been recorded from four of the Yilgarn BIFs: 
Koolyanobbing (KOO), Windarling (WIN), Mt Jackson (MTJ) 
and the Die Hardy Range (DH) (pers. comm. R. Teale) [33]. This 
study extended the known distribution to include the Helena 
Aurora (HA) Range. Four more BIFs in the region (Mt Dimer, Mt 
Correll, Hunt Range and Mt Manning) were also searched but no 
specimens were found. Specimens from the Windarling and 
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Koolyanobbing Ranges, preserved in 100% ethanol were available 
from the Western Australian Museum's Department of Terrestrial 
Zoology collection (Table SI). For the remaining BIFs (DH, MTJ 
and HA) sampling was spread across the geographic extent of each 
BIF as far as was practical to avoid sampling closely related 
individuals. For some BIFs this was not possible as millipedes were 
either not found and/ or there was no suitable habitat. The same 
search procedure was applied to all BIFs, so we believe it is 
reasonable to assume that the specimens collected reflect the 
genetic diversity contained within each site. The five BIFs vary in 
size, topographic complexity and orientation. Koolyanobbing is 
large, approximately 12 km in length and reaching heights of 
500 m. This BIF is topographically complex and contains many 
sites with habitat suitable for millipedes. Helena Aurora is also 
large (~ 12 km) and elevated (600-700 m height), but is comprised 
of more rounded hills with gentle slopes. Millipedes here were 
associated with drainage lines coming off the peaks. The Die 
Hardy Range is large and X-shaped, a topology providing many 
shaded slopes that supported habitat suitable for millipedes. The 
two smaller Ranges Windarling and Mt Jackson, have a narrow, 
east to west orientation and are therefore subject to an abundance 
of northern sunlight. These smaller Ranges did not possess as 
many water-gaining sites and millipedes were scarcer. Invertebrate 
surveys conducted on the intervening 'flats' between the Yilgarn 
BIF have not recorded A. bamfordi (pers. comm. R. Teale). We also 
searched any water-gaining or densely vegetated sites between 
BIFs and although these recovered other millipedes from the genus 
Antichiropus, A. bamfordi was not found. Thus, to the best of our 
knowledge, A. bamfordi populations are restricted to BIF habitat 
and are not presently connected by populations in intervening 
habitat. The field study did not involve any endangered or 
protected species and none of the field sites from which we 
sampled (i.e. excluding KOO and WIN) required access permis- 
sion, with the exception of Mt Jackson (permission required from 
Cliffs Natural Resources) and Mt Manning (Regulation 4 permit 
required from the Western Australian Department of Parks and 
Wildlife (DPaW)). A license to take fauna for scientific purposes 
was obtained prior to sampling from DPaW. GPS coordinates of 
the sampling locations are provided in Table S 1 . 

In total, 77 individuals from five BIFs were sequenced at two 
mitochondrial regions (77 at 16 S, 75 at COl) and genotyped at 
1 1 nuclear microsatellite loci. DNA was extracted from the legs or 
from the crushed heads of juvenile specimens stored in 100% 
ethanol using a Qiagen DNeasy Kit, or a salt-based DNA 
extraction protocol [40]. Elutions were performed in 100 (il or 
50 ul of TE buffer respectively. Legs and head tissue were 
homogenized in the lysis buffer prior to addition of proteinase K 
using mini micropestles (Geneworks). 

Mitochondrial DNA amplification and sequencing 

We sequenced two partial regions of the mitochondria, the 
Cytochrome Oxidase Subunit 1 (COl) and the large Ribosomal 
Subunit 16 S. For amplification of COl we used the Folmer 
primers LCO1490 (5 ' GGTC AAC AAATC ATAAAGATAT- 
TGG3') and HC02198 (5'TAAACTTCAGGGTGACCAAA- 
AAATCA3') [41] and for 16 S, the 'universal' primers 16Sar 
(5 ' CGCCTGTTTTTC AAAAACAT3 ') and 16Sbr (5'CCGGTT- 
TGAACTCAGATCATGT3') [42]. We also amplified two 
nuclear gene regions, the Internal Transcribed Spacer 1 (ITS-1) 
[43] and the 18S rRNA gene [44]. ITS-1 sequenced products 
showed evidence of site heterogeneity and were excluded from use, 
and 18 S sequences were monomorphic. Amplifications were 
performed in 25 ul reactions containing 1 x PCR buffer, 3 mM 
MgCl 2 , 200 uM dNTP's, 0.4 uM of each forward and reverse 



primer, 0.25 u.g of Bovine Serum Albumin (Fisher Biotech), 0.3 U 
of Platinum Taq DNA Polymerase (Invitrogen) and approximately 
20 ng/ul DNA. The thermal cycler profile for both COl and 16 S 
involved an initial denaturation at 94°C for 3 min, followed by 35 
cycles of 94°C for 30 sec, 46°C for 30 sec and 72°C for 30 sec 
followed by a final extension at 72°C for 2 min. PCR products 
were purified using an UltraClean PCR Cleanup kit (Mo-Bio 
technologies, Geneworks) and the cleaned product was quantified 
using a Nanodrop (Nanodrop Technologies). Sequencing reactions 
(1/8) were carried out in 10 (il volumes containing 1 of Big Dye 
Terminator, 1.5 ul of 5 x sequencing buffer, 1 ul of 3.2 pmol 
primer and approximately 40 ng of purified PCR product. 
Thermal cycling conditions for sequencing reactions involved an 
initial denaturation at 96°C for 2 min, followed by 25 cycles of 
96°C for 10 sec, 46°C for 5 sec and 60°C for 4 min. Sequenced 
products were cleaned prior to electrophoresis using a standard 
ethanol precipitation method (Applied Biosystems). Products were 
sequenced on an Applied Biosystems 3730 capillary sequencer 
(Murdoch University). 

Mitochondrial DNA analysis 

Sequence data were aligned using Clustal W in the program 
BioEdit [26] and finalised by eye with all mutations checked 
against raw sequence data. We used AMOVA, implemented in 
Arlequin v. 3.5 [45] to determine the structuring of mitochondrial 
DNA diversity amongst BIF populations using the number of 
pairwise base-pair differences and tested for significance with 
10000 permutations. Bayesian analysis was used to visualise 
genetic relationships and to estimate coalescence times of the 
populations in BEAST v. 1.7.5 [46]. Data were partitioned into 
the COl and 16 S datasets, enforcing monophyly for all A. 
bamfordi samples, with a SRD06 codon-partitioned model of 
sequence evolution applied to the COl dataset and an HKY+G 
model to the 16 S partition. Due to an absence of fossil data and 
geological calibration points, estimates of population coalescence 
times were generated by applying the average arthropod mutation 
rate of 2.3% per million years to the partitioned data set [47]. We 
applied a strict clock to both partitions and used a Coalescent- 
Constant Size tree prior for four independent runs of 40 million 
generations, sampling every 1000 generations. Log files were 
combined in Log Combiner and the first 25% discarded as burnin. 
Visual inspection of the traces and high ESS values (>200) 
confirmed stationarity had been reached in Tracer (BEAST 
package). Tree files were combined in LogCombiner and the final 
file annotated in TreeAnnotater (BEAST package). 

To distinguish between incomplete lineage sorting and migra- 
tion for the two populations that were paraphyletic (see results), we 
performed nested model likelihood ratio tests of various isolation 
versus migration models in 'L mode' in the program IMa2 [48]. 
IMa2 uses Markov chain Monte Carlo simulations of gene 
genealogies to estimate divergence times of populations and 
effective population sizes of extant and ancestral populations [49] . 
To assess whether migration or drift was responsible for the 
paraphyly observed in the two populations we performed a run on 
the two populations in 'M mode' using the following command 
line priors: q200, tlO, ml, b 100000, 11000000, hfg, hn6, ha0.96, 
hb0.9. Convergence was assessed by monitoring ESS values (all 
>20 000), chain swap rates (>80%) and running the same 
command line four times independendy using different seed values 
to ensure consistency of results. The output was then used to test 
the likelihood of 24 different migration/isolation models oudined 
in Hey and Nielsen [50] using the following command line: q200, 
tlO, ml, rO, c2 b 100000 11000000. Likelihood ratio tests were used 
to reject models that were significantly worse than the full model, 
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and the remaining models were compared using the AIG criterion 
(AIC = -2(log(P)+2K)) to identify the model of best fit, where K is 
the number of parameters in the statistical model and L is the 
maximised value of the likelihood function [51]. 

To assess the historical stability of BIF population sizes we 
performed tests of neutrality that can also serve as a measure of 
demographic expansion or contraction. Estimates of Tajima's D 
were calculated in Arlequin. Diversity statistics including haplo- 
type diversity, nucleotide diversity (Pi), and nucleotide divergence 
Jukes and Cantor) between populations were calculated in DnaSP 
v. 5.10 [52]. 

To test for isolation by distance within BIF populations using 
the combined sequence data, we conducted mantel tests in R using 
the Vegan package [53]. Permutations were performed 10000 
times on matrices of the number of pairwise sequence differences 
and geographic distance. The R statistic is based on Pearson's 
product-moment correlation coefficient. We plotted the number of 
pairwise differences on geographic distance and established a 
linear relationship meaning no transformations of the data were 
required. 

Nuclear microsatellite genotyping and data analysis 

Genotyping of the 1 1 microsatellite markers was carried out 
according to Nistelberger et al. [54]. We checked all pairwise 
combinations of loci for linkage disequilibrium by performing 
exact tests in Genepop v. 1.2 [55] and applying a sequential 
Bonferroni correction to determine significance [56]. There was 
no significant association among loci following sequential 
Bonferroni correction. Traditional F statistics were of limited use 
in our study due to the likelihood that our BIF 'populations' do not 
represent true panmictic populations given the distances over 
which individuals were sampled and the limited dispersal 
capabilities of A. bamfordi. This was confirmed with all populations 
showing significant deviation from Hardy- Weinberg equilibrium, 
with the exception of Windarling, where individuals had been 
sampled in close proximity to one another. Rarefied allelic 
richness provided a better estimator of genetic diversity, given the 
sampling, and was calculated in HP-Rare v 1.0 [57]. The 
partitioning of genetic diversity amongst BIF populations was 
assessed using an Analysis of Molecular Variance (AMOVA) 
implemented in Arlequin with significance determined with 1 0000 
permutations. An analysis of genetic structure across the five BIF 
populations was assessed using a Bayesian clustering algorithm 
implemented in TESS ver. 2.3 [58]. This program assumes 
proximate interactions between individuals with spatial informa- 
tion prescribed at the individual level. We ran TESS with no 
admixture, owing to the limited dispersal capabilities of millipedes 
and performed 100 iterations of K values between 1 and 7, with a 
burnin length of 20000, 100000 sweeps and a spatial interaction 
parameter of 0.8. Lowest DIC values for each K value were 
plotted in order to determine the optimal K value. The results 
were then used in the downstream program CLUMPP [59] which 
aligns cluster assignment across the replicate analyses and ensures 
the correct assignment of K. Structure output graphics were 
created using DISTRUCT [60]. 

To further reconstruct the phylogenetic relationships of 
populations using microsatellite data we generated Cavalli-Sforza 
chord distances (D c ) [61] between populations and determined 
support using 100 bootstrap replications in the program Popula- 
tions 1.2.30 [62]. This geometric distance measure does not make 
any biological assumptions about the data and can be more useful 
for phylogenetic inference than estimates that use stepwise 
mutation models when small population sizes are involved 
[61,63]. Distances were imported into TreeView for visualisation 



[64] . To assess isolation by distance within BIF populations using 
microsatellite data we correlated pairwise population matrices of 
genetic distance and geographic distance and tested for signifi- 
cance using a Mantel test with 9999 permutations in the program 
IBDWS v 3.23 [65]. 

Results 

Genetic variation 

Mitochondrial DNA data. Of the two mitochondrial regions 
amplified, COl was the most variable, with 45 haplo types 
identified from the 75 individuals sequenced (two WIN samples 
did not amplify at COl). The variation was composed of 54 
transitions, eight transversions and two multistate substitutions. 
The 16 S region yielded 24 haplotypes from the 77 individuals 
sequenced and included 21 transitions, one transversion and one 
single base pair indel. Combining the two regions resulted in a 
total, aligned sequence length of 1 1 14 bp and 48 haplotypes. The 
most prevalent haplotype (Hap 47 WIN) had a frequency of 23%, 
followed by Hap 35 (MTJ) and Hap 39 (HA) with frequencies of 
8%, and Hap 33 (MTJ) and Hap 40 (HA) with 6%. The remaining 
haplotypes were only found in one or two individuals. Measures of 
genetic diversity were consistently lowest in the Windarling 
population and highest in Koolyanobbing. Haplotype diversity 
ranged from 0.295 (WIN) to 0.986 (KOO) with an average of 
0.794. Nucleotide diversity ranged from 0.06% (WIN) to 8.5% 
(KOO) with an average of 4%. No populations departed from the 
neutral model of evolution except for Windarling (D = —1.775) 
(Table 1). Measures of diversity in Windarling are likely to be 
biased due to restricted sampling. Millipedes had been previously 
collected from two, closely located sites on this BIF (<260 m 
apart), and although noted to be rare along the length of the 
Range (pers. comm. R. Teale), further sampling was not possible 
due to restricted access associated with mining activity. 

There was a large and highly significant degree of genetic 
structuring of the mitochondrial genetic diversity, with most 
variation (AMOVA) occurring among BIF populations (71.56%) 
rather than within (28.4%). The greatest genetic divergence 
occurred between the Mt Jackson and Helena Aurora populations 
(approx. 40 km apart) and the least between the two most spatially 
isolated populations, Koolyanobbing and Die Hardy (approx. 
100 km apart) (Table 2, Fig. 1). 

Nuclear microsatellite data. Microsatellite genetic diversity 
was low, both within and amongst BIF populations. Rarefied 
allelic richness ranged from 1.89 (HA) to 3.28 (KOO) with an 
average of 2.47. There was a highly significant degree of 
structuring of the microsatellite diversity with the majority of 
variation (68. 1 1 %) maintained within populations, compared to 
(31.89%) among populations. 

Phylogeographic structure and population genetic 
structure 

There was evidence of strong phylogeographic structure within 
A. bamfordi, with minimal sharing of COl and 16 S haplotypes 
amongst individuals within BIF populations and no sharing of 
haplotypes between populations. Bayesian analysis identified five 
genetic clades in which populations were reciprocally monophy- 
letic, with the exception of Koolyanobbing and Helena Aurora for 
which some samples were paraphyletic (Fig. 2). This was due to a 
subset of six individuals from the north of the Koolyanobbing 
Range being grouped with samples from Helena Aurora (Fig. 1). 
Support for the five genetic clades was high (posterior probabilities 
^0.85) although the relationships between the clades was poorly 
supported. Dating estimates showed coalescence during the 
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Figure 1. Schematic representation of the five Atelomastix 
bamfordi Banded Iron Formation populations, coloured ac- 
cording to the five genetic clades identified with mtDNA 
analysis. Black circles indicate sample sites. Study area shown in inset. 
doi:10.1371/joumal.pone.0093038.g001 

Pleistocene, between 0.729 Ma [0.532-0.952] and 1.078 Ma 
[0.792-1.433], although only the date of 1.078 for the divergence 
of the Mt Jackson population has high support (Fig. 2). 

Isolation and migration analysis in the program IMa2 was used 
to determine whether the paraphyly between Koolyanobbing and 
Helena Aurora was due to incomplete lineage sorting (isolation) or 
gene flow (migration). The three models with the highest AIC 
support all indicated there was no historical migration between the 
two populations (Table S2). 

Mantel tests of pairwise sequence differences against geographic 
distance indicated significant isolation by distance within the 
Koolyanobbing (r = 0.740 p = 0.0001), Mt Jackson (r = 0.225 
p = 0.0375) and Helena Aurora (r = 0.399 p = 0.0025) populations. 

Population genetic structure 

Bayesian cluster analysis identified four genetic clusters in the 
microsatellite data. The first corresponded to the Die Hardy and 
Mt Jackson populations, and the remaining clusters to each of the 
three other populations (HA, WIN, KOO) (Fig. 3). Further 
investigation of structure in the first cluster revealed two more 
groupings corresponding to the separate BIF populations (DH, 
MTJ). There was some evidence of admixture in the Mt Jackson 
and Koolyanobbing populations (Fig. 3). Phylogenetic analysis of 
relationships between populations using Cavalli-Sforza chord 
distances showed distinct genetic structure amongst populations 
using microsatellite data, but with some clustering of Koolyanob- 
bing and Helena Aurora (97% bootstrap support) and Die Hardy 
and Mt Jackson (62%) (Fig. 4). Two of the five populations, Helena 



Table 2. Pairwise population nucleotide divergence based on 
Jukes and Cantor of five Atelomastix bamfordi populations 
using mtDNA data. 





DH 


HA 


MTJ 


WIN 


KOO 


DH 


0 










HA 


0.016 


0 








MTJ 


0.017 


0.022 


0 






WIN 


0.013 


0.018 


0.019 


0 




KOO 


0.013 


0.014 


0.018 


0.017 


0 



*P value<0.005. 

doi:1 0.1 371 /journal.pone.0093038.t002 
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Mt Jackson 

Koolyanobbing 
Koolyanobbing 
Helena Aurora 



Die Hardy 



Windarling 



Figure 2. Bayesian cladogram of 75 Atelomastix bamfordi individuals based on combined COI and 16 S mtDNA data. Clades are 
designated separate colours as per Fig. 1 . Numbers on terminals indicate haplotype number. The * denotes major nodes with high support (posterior 
probability>0.85) and the numbered node indicates divergence estimate in millions of years with confidence intervals in brackets. NB the A. 
nigrescens branch has been truncated (//) for viewing purposes. 
doi:1 0.1 371 /journal.pone.0093038.g002 





DH MTJ HA WIN KOO 

Figure 3. Bayesian population assignment of Atelomastix 
bamfordi individuals generated in TESS for the optimal 
grouping of K = 4. Each vertical line represents an individual and 
the colour displays the probability of belonging to each genetic cluster. 
doi:1 0.1 371 /journal.pone.0093038.g003 



KOO 




WIN 



DH 

Figure 4. Unrooted neighbour-joining tree of Cavalli-Sforza 
chord distances between the five BIF populations of Atelomas- 
tix bamfordi. Distances are based on 1 1 microsatellite loci with 
bootstrap support shown at nodes. 
doi:10.1371/journal.pone.0093038.g004 
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Aurora (r = 0.30 p = 0.004) and Die Hardy (r = 0.32 p = 0.009), 
showed weak but significant patterns of isolation-by-distance at 
microsatellite loci following Mantel tests. 

Discussion 

Genetic differentiation of BIF populations 

Spatial patterns. Phylogeographic analysis of A. bamfordi has 
revealed strong genetic differentiation between BIF populations, 
with mtDNA genetic lineages corresponding to each of the five 
BIF populations except for six individuals from Koolyanobbing 
that were paraphyletic with the Helena Aurora Range due to 
incomplete lineage sorting. Microsatellite data were broadly 
consistent with this result, identifying four genetic clusters, the 
first of which was further broken down into individual populations. 

The lack of correlation between genetic divergence and 
geographic distance, with the two most spatially isolated popula- 
tions showing the least genetic divergence, suggests there has been 
littie, to no connectivity of BIF populations following their 
separation. This result was emphasized in the Bayesian analysis 
which showed high support for genetic clades but poor support for 
the relationships between clades (Fig. 2). It is also consistent with 
estimates of sequence neutrality (D) that indicated a lack of 
demographic expansion and contraction events, with the excep- 
tion of Windarling, where the result can be explained by close 
sample proximity. This pattern of differentiation has been 
observed in other taxa with restricted dispersal capabilities where 
populations have become isolated as a response to climate change 
and habitat loss. Like A. bamfordi, these studies emphasise the role 
of genetic drift in driving population differentiation. For example, 
mitochondrial and nuclear sequence data revealed extraordinarily 
strong phylogeographic structure in populations of trapdoor 
spiders (Moggridgea tingle) in south-western Australia, often situated 
just kilometres apart. Molecular dating suggested these popula- 
tions became isolated in the mesic valleys and uplands of the 
Porongorup and Stirling Ranges following the onset of aridity in 
the late Miocene [66] . Like A. bamfordi, reciprocally monophyletic 
lineages across three geographic locations highlighted the repro- 
ductive isolation of M. tingle populations but the pattern 
emphasized a more ancient division, with monophyly observed 
in both mtDNA and the more slowly evolving nuclear genome 
[66]. Studies of genetic structure in two frog species in south- 
western Australia, Geocrinia lutea and G. rosea, also described 
patterns of extremely low dispersal and long-term isolation 
between habitat 'islands', with divergence attributed to isolation 
via range changes resulting from climate change rather than 
geographical barriers [67]. Another striking example of genetic 
divergence across terrestrial islands occurred in the sky islands of 
south-eastern Arizona - forested mountain Ranges, separated by 
low-lying grassland and desert. Here, the jumping spider 
Habronattus pugillis occupies woodland habitat at altitudes above 
1400 m but is absent from the lowlands. Using mtDNA, Masta [7] 
showed that like A. bamfordi, smaller populations of H. pugillis 
formed monophyletic clades whereas larger Ranges showed 
evidence of paraphyly. Like A. bamfordi, this was ascribed to 
incomplete lineage sorting rather than resulting from migration 
[7] . The isolation of populations was attributed to climate-induced 
contraction of woodland populations onto disjunct mountain 
Ranges. 

Temporal patterns. Molecular dating with the arthropod 
clock rate has placed the separation of A. bamfordi clades at 0.7 to 
1.1 million years ago during the Pleistocene. The timing coincides 
with known increases in aridity in the transitional rainfall zone 
(TRZ) due to climatic fluctuations associated with Milankovitch 



cycles during this period [68]. In particular there was a major arid 
shift in southern Australia corresponding to the middle Pleistocene 
transition (1250 ka-700 ka), a period of increased severity of 
glaciation, associated with greater ice volumes, lowered temper- 
atures and increased aridity [69] . Although the TRZ remained un- 
glaciated, this period of intense climatic instability is expected to 
have caused the fragmentation of species distributions due to 
changes in suitable habitat within this region [70]. Although the 
divergence estimates presented here must be treated with caution 
given the use of a generalised arthropod mutation rate, the errors 
associated with this method are lessened given we are working at 
shallow phylogenetic levels [71,72]. 

BIF population dynamics - historical and contemporary 
patterns 

Historical patterns. The process of genetic drift acting 
independendy on each BIF population is reflected in different 
patterns of genetic divergence and diversity seen in the mtDNA 
data. The smaller populations, Mt Jackson and Windarling, and 
the larger population Helena Aurora, which had less suitable 
millipede habitat, showed lower genetic diversity and higher levels 
of genetic divergence than the larger populations (KOO and DH). 
These three BIFs have probably harboured historically small 
effective population sizes, owing to the limited availability of moist 
sites, increasing their susceptibility to forces of genetic drift driving 
divergence [73]. They may also be influenced by different selective 
forces associated with the more arid conditions present on these 
BIFs. On the other hand the two large BIF populations, 
characterised by a greater availability of suitable habitat (KOO 
and DH), retained genetic diversity and had less genetic 
divergence, consistent with historically large population sizes 
[73]. The paraphyly of individuals from Koolyanobbing with those 
of Helena Aurora further supports the likelihood of a large 
historical population size at Koolyanobbing, as lineage sorting is 
still ongoing at this BIF. Over time this signature is likely to erode, 
leaving the five BIFs reciprocally monophyletic. 

Contemporary patterns. Rapidly evolving microsatellite 
markers can provide a more contemporary view of BIF population 
dynamics. The information we gained from these markers was 
limited by a lack of Hardy- Weinberg equilibrium within our BIF 
populations yet we still identified patterns in genetic diversity that 
appear to support the impact of genetic drift on populations: lower 
genetic diversity in BIF's likely to have had small effective 
population sizes (WIN, MTJ and HA). Signals of isolation by 
distance (IBD) were seen in some of the BIF populations (HA, 
DH), but these were weak and likely to be heavily influenced by 
limited sample sizes. The small pairwise distances between samples 
at Windarling explains the lack of IBD seen in this BIF, but we 
were surprised to see evidence of IBD in the Helena Aurora and 
Mt Jackson individuals. Here, we expected limited gene flow to 
result in a pattern of 'islands within islands', where significant 
isolation has resulted in genetically divergent populations within 
BIFs that do not correlate with geographic distance. An increased 
number of samples from more locations would be required to 
determine whether IBD is really occurring at these spatial scales. 

Conservation implications 

The evolutionary significance of geographically isolated genetic 
lineages has long been recognised and conservation biologists have 
attempted to formalise this by delineating conservation units 
within a species such as Evolutionary Significant Units (ESUs) and 
Management Units (MUs) [74]. Conservation units have been 
defined at various levels of genetic differentiation within species 
and are likely to experience limited gene flow or even total 
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reproductive isolation. For example, ESUs have been defined as a 
population or a group of populations that warrant separate 
management for conservation because of high genetic and 
ecological distinctiveness [75], with Mortiz [20] suggesting that 
ESUs might be characterised by a pattern of reciprocal monophyly 
amongst populations. Applying this criterion to our mtDNA data, 
results in the recognition of three BIF populations (DH, MTJ, 
WIN) as ESUs. However, Paetkau [76] and Crandall et al. [77] 
highlight the problems inherent in ignoring paraphyletic lineages, 
such as Koolyanobbing and Helena Aurora, which are often 
ancestral and maintain high levels of diversity and divergence 
within themselves. In our case, the paraphyly of Koolyanobbing 
with Helena Aurora is a result of a slower progression towards 
monophyly due to ongoing lineage sorting, probably owing to 
large population sizes. Given the clear evidence for restricted gene 
flow between these two populations, both could readily be 
recognised as MUs. Therefore we consider that all five BIF 
populations warrant designation as conservation units. 

Although it is difficult to predict future climate change at such 
fine scales, broad-scale models for this region over the 21 st century 
indicate increasing temperature, decreasing rainfall and a greater 
incidence of fire [78]. If these predictions are correct, populations 
are unlikely to re-connect via suitable intervening habitat in the 
future, further consolidating their isolation to the remaining moist 
pockets provided by BIF, where they will continue to be influenced 
by forces of genetic drift and local selection. Given their 
dependence on moisture, conservation efforts should focus on 
protecting areas of dense leaf litter, large outcrops and shaded 
gullies on each BIF in order to preserve the five genetic lineages 
identified here. To achieve this, preservation of the height and 
topographical complexity of BIFs is essential, in order to maintain 
drainage lines and water-collecting sites. Conservation of these 
areas is also likely to support populations of other, relict fauna and 
flora that are adapted to moist environments. 

Conclusion 

This is the first phylogeographic study conducted on a species 
inhabiting multiple Banded Iron Formations in this landscape and 
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